<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of VMT_ProcessTransectsV3_new</title>
  <meta name="keywords" content="VMT_ProcessTransectsV3_new">
  <meta name="description" content="Not used by current version">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2005 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../m2html.css">
  <script type="text/javascript">
    if (top.frames.length == 0) { top.location = "../../index.html"; };
  </script>
</head>
<body>
<a name="_top"></a>
<!-- ../menu.html VMT_4.0_dev --><!-- menu.html depreciated -->
<h1>VMT_ProcessTransectsV3_new
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>Not used by current version</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>function [A,V] = VMT_ProcessTransectsV3_new(z,A,setends) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre class="comment"> Not used by current version
This routine acts as a driver program to process multiple transects at a
single cross-section for velocity mapping.</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="VMT_GridData2MeanXS.html" class="code" title="function [A,V] = VMT_GridData2MeanXS(z,A,V)">VMT_GridData2MeanXS</a>	Generates a uniformly spaced grid for the mean cross section and</li></ul>
This function is called by:
<ul style="list-style-image:url(../../matlabicon.gif)">
</ul>
<!-- crossreference -->

<h2><a name="_subfunctions"></a>SUBFUNCTIONS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="#_sub1" class="code">function [A,V] = VMT_MapEns2MeanXSV2(z,A,setends)</a></li><li><a href="#_sub2" class="code">function [A,V] = VMT_GridData2MeanXS(z,A,V)</a></li><li><a href="#_sub3" class="code">function mypostcallback_zoom(obj,evd)</a></li><li><a href="#_sub4" class="code">function mypostcallback_pan(obj,evd)</a></li><li><a href="#_sub5" class="code">function [A,V] = VMT_CompMeanXS_old(z,A,V)</a></li><li><a href="#_sub6" class="code">function [A,V] = VMT_CompMeanXS_UVW(z,A,V)</a></li><li><a href="#_sub7" class="code">function [A,V] = VMT_CompMeanXS_PriSec(z,A,V)</a></li><li><a href="#_sub8" class="code">function [V] = VMT_RozovskiiV2(V,A)</a></li></ul>

<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [A,V] = VMT_ProcessTransectsV3_new(z,A,setends)</a>
0002 <span class="comment">% Not used by current version</span>
0003 <span class="comment">%This routine acts as a driver program to process multiple transects at a</span>
0004 <span class="comment">%single cross-section for velocity mapping.</span>
0005 
0006 <span class="comment">%V2 adds the cpability to set the endpoints of the mean cross section</span>
0007 
0008 <span class="comment">%V3 adds the Rozovskii computations for secondary flow. 8/31/09</span>
0009 
0010 <span class="comment">%Among other things, it:</span>
0011 
0012 <span class="comment">% Determines the best fit mean cross-section line from multiple transects</span>
0013 <span class="comment">% Map ensembles to mean c-s line</span>
0014 <span class="comment">% Determine uniform mean c-s grid for vector interpolating</span>
0015 <span class="comment">% Determine location of mapped ensemble points for interpolating</span>
0016 <span class="comment">% Interpolate individual transects onto uniform mean c-s grid</span>
0017 <span class="comment">% Average mapped mean cross-sections from individual transects together</span>
0018 <span class="comment">% Rotate velocities into u, v, and w components</span>
0019 
0020 
0021 <span class="comment">%(adapted from code by J. Czuba)</span>
0022 
0023 <span class="comment">%P.R. Jackson, USGS, 12-9-08</span>
0024 
0025 disp(<span class="string">'Processing Data...'</span>)
0026 <span class="comment">% warning off</span>
0027 <span class="comment">%% Map ensembles to mean cross-section</span>
0028 
0029 [A,V] = <a href="#_sub1" class="code" title="subfunction [A,V] = VMT_MapEns2MeanXSV2(z,A,setends)">VMT_MapEns2MeanXSV2</a>(z,A,setends);
0030 
0031 <span class="comment">%% Grid the measured data along the mean cross-section</span>
0032 <span class="comment">%[A,V] = VMT_GridData2MeanXS(z,A,V);</span>
0033 [A,V] = <a href="VMT_GridData2MeanXS.html" class="code" title="function [A,V] = VMT_GridData2MeanXS(z,A,V)">VMT_GridData2MeanXS</a>(z,A,V);
0034 
0035 <span class="comment">%% Computes the mean data for the mean cross-section</span>
0036 <span class="comment">%[A,V] = VMT_CompMeanXS(z,A,V);</span>
0037 [A,V] = <a href="#_sub5" class="code" title="subfunction [A,V] = VMT_CompMeanXS_old(z,A,V)">VMT_CompMeanXS_old</a>(z,A,V);
0038 
0039 <span class="comment">%% Decompose the velocities into u, v, and w components</span>
0040 [A,V] = <a href="#_sub6" class="code" title="subfunction [A,V] = VMT_CompMeanXS_UVW(z,A,V)">VMT_CompMeanXS_UVW</a>(z,A,V);
0041 
0042 <span class="comment">%% Decompose the velocities into primary and secondary components</span>
0043 [A,V] = <a href="#_sub7" class="code" title="subfunction [A,V] = VMT_CompMeanXS_PriSec(z,A,V)">VMT_CompMeanXS_PriSec</a>(z,A,V);
0044 
0045 <span class="comment">%% Perform the Rozovskii computations</span>
0046 V = <a href="#_sub8" class="code" title="subfunction [V] = VMT_RozovskiiV2(V,A)">VMT_RozovskiiV2</a>(V,A);
0047 
0048 <span class="comment">%%</span>
0049 <span class="comment">%figure(4); clf</span>
0050 <span class="comment">%plot3(V.mcsX,V.mcsY,V.mcsDepth(1))</span>
0051 disp(<span class="string">'Processing Completed'</span>)
0052 
0053 <span class="comment">%% Notes:</span>
0054 
0055 <span class="comment">%1. I removed scripts at the end of the original code that computes</span>
0056 <span class="comment">%standard deviations (12-9-08)</span>
0057 
0058 <span class="comment">% else</span>
0059 <span class="comment">%</span>
0060 <span class="comment">% disp('Processing Data...')</span>
0061 <span class="comment">% %% Map ensembles to mean cross-section</span>
0062 <span class="comment">% hf = VMT_CreatePlotDisplay('shiptracks');</span>
0063 <span class="comment">%</span>
0064 <span class="comment">% [A,V] = VMT_MapEns2MeanXSV2_PDF(hf,z,A,setends);</span>
0065 <span class="comment">%</span>
0066 <span class="comment">% %% Grid the measured data along the mean cross-section</span>
0067 <span class="comment">% %[A,V] = VMT_GridData2MeanXS(z,A,V);</span>
0068 <span class="comment">% [A,V] = VMT_GridData2MeanXS_PDF(hf,z,A,V);</span>
0069 <span class="comment">%</span>
0070 <span class="comment">% return</span>
0071 <span class="comment">% %% Computes the mean data for the mean cross-section</span>
0072 <span class="comment">% %[A,V] = VMT_CompMeanXS(z,A,V);</span>
0073 <span class="comment">% [A,V] = VMT_CompMeanXS_old(z,A,V);</span>
0074 <span class="comment">%</span>
0075 <span class="comment">% %% Decompose the velocities into u, v, and w components</span>
0076 <span class="comment">% [A,V] = VMT_CompMeanXS_UVW(z,A,V);</span>
0077 <span class="comment">%</span>
0078 <span class="comment">% %% Decompose the velocities into primary and secondary components</span>
0079 <span class="comment">% [A,V] = VMT_CompMeanXS_PriSec(z,A,V);</span>
0080 <span class="comment">%</span>
0081 <span class="comment">% %% Perform the Rozovskii computations</span>
0082 <span class="comment">% V = VMT_RozovskiiV2(V,A);</span>
0083 <span class="comment">%</span>
0084 <span class="comment">% %%</span>
0085 <span class="comment">% %figure(4); clf</span>
0086 <span class="comment">% %plot3(V.mcsX,V.mcsY,V.mcsDepth(1))</span>
0087 <span class="comment">% disp('Processing Completed')</span>
0088 <span class="comment">%</span>
0089 <span class="comment">% %% Notes:</span>
0090 <span class="comment">%</span>
0091 <span class="comment">% %1. I removed scripts at the end of the original code that computes</span>
0092 <span class="comment">% %standard deviations (12-9-08)</span>
0093 <span class="comment">%</span>
0094 <span class="comment">% end</span>
0095 
0096 <span class="comment">%==========================================================================</span>
0097 <a name="_sub1" href="#_subfunctions" class="code">function [A,V] = VMT_MapEns2MeanXSV2(z,A,setends)</a>
0098 
0099 <span class="comment">%This routine fits multiple transects at a single location with a single</span>
0100 <span class="comment">%line and maps individual ensembles to this line. Inputs are number of files (z) and data matrix (Z)(see ReadFiles.m).</span>
0101 <span class="comment">%Output is the updated data matrix with new mapped variables.</span>
0102 
0103 <span class="comment">%V2 adds the capability to set the endpoints of the mean cross section</span>
0104 
0105 <span class="comment">%(adapted from code by J. Czuba)</span>
0106 
0107 <span class="comment">%P.R. Jackson, USGS, 12-9-08</span>
0108 
0109 <span class="comment">%% Determine the best fit mean cross-section line from multiple transects</span>
0110 <span class="comment">% initialize vectors for concatenation</span>
0111 
0112 x = [];
0113 y = [];
0114 <span class="comment">% figure(1); clf</span>
0115 <span class="comment">% set(gca,'DataAspectRatio',[1 1 1],'PlotBoxAspectRatio',[1 1 1])</span>
0116 mfd = zeros(1,z);
0117 <span class="keyword">for</span> zi = 1 : z
0118     
0119     <span class="comment">% concatenate coords into a single column vector for regression</span>
0120     x = cat(1,x,A(zi).Comp.xUTM);
0121     y = cat(1,y,A(zi).Comp.yUTM);
0122 
0123 <span class="comment">%     figure(1); hold on</span>
0124 <span class="comment">%     %plot(A(zi).Comp.xUTM,A(zi).Comp.yUTM,'r'); hold on</span>
0125 <span class="comment">%     plot(A(zi).Comp.xUTMraw,A(zi).Comp.yUTMraw,'b'); hold on</span>
0126     
0127     <span class="comment">%Find the mean flow direction for each transect</span>
0128     mfd(zi) = nanmean(nanmean(A(zi).Clean.vDir,1)); 
0129 <span class="keyword">end</span>
0130 V.mfd = nanmean(mfd); <span class="comment">% Mean flow direction for all the transects</span>
0131 
0132 <span class="keyword">if</span> setends  <span class="comment">%Gets a user text file with fixed cross section end points</span>
0133     
0134     defaultpath = <span class="string">'C:\'</span>;
0135     defaultpath = pwd;
0136     endspath = [];
0137     <span class="keyword">if</span> exist(<span class="string">'\VMT\LastDir.mat'</span>) == 2
0138         load(<span class="string">'VMT\LastDir.mat'</span>);
0139         <span class="keyword">if</span> exist(endspath) == 7
0140             [file,endspath] = uigetfile({<span class="string">'*.txt;*.csv'</span>,<span class="string">'All Text Files'</span>; <span class="string">'*.*'</span>,<span class="string">'All Files'</span>},<span class="string">'Select Endpoint Text File'</span>,endspath);       
0141         <span class="keyword">else</span>
0142             [file,endspath] = uigetfile({<span class="string">'*.txt;*.csv'</span>,<span class="string">'All Text Files'</span>; <span class="string">'*.*'</span>,<span class="string">'All Files'</span>},<span class="string">'Select Endpoint Text File'</span>,defaultpath);
0143         <span class="keyword">end</span>
0144     <span class="keyword">else</span>
0145         [file,endspath] = uigetfile({<span class="string">'*.txt;*.csv'</span>,<span class="string">'All Text Files'</span>; <span class="string">'*.*'</span>,<span class="string">'All Files'</span>},<span class="string">'Select Endpoint Text File'</span>,defaultpath);
0146     <span class="keyword">end</span>
0147     
0148     <span class="keyword">if</span> ischar(file)
0149         infile = [endspath file];
0150         <span class="comment">%[file,path] = uigetfile({'*.txt;*.csv','All Text Files'; '*.*','All Files'},'Select Endpoint Text File');</span>
0151         <span class="comment">%infile = [path file];</span>
0152         disp(<span class="string">'Loading Endpoint File...'</span> );
0153         disp(infile);
0154         data = dlmread(infile);
0155         x = data(:,1);
0156         y = data(:,2);
0157         figure(1); hold on
0158         plot(x,y,<span class="string">'go'</span>,<span class="string">'MarkerSize'</span>,10); hold on
0159     <span class="keyword">end</span>
0160 <span class="keyword">end</span> 
0161 
0162 <span class="comment">% find the equation of the best fit line</span>
0163 xrng = max(x) - min(x);
0164 yrng = max(y) - min(y);
0165 
0166 <span class="keyword">if</span> xrng &gt;= yrng <span class="comment">%Fit based on coordinate with larger range of values (original fitting has issues with N-S lines because of repeated X values), PRJ 12-12-08</span>
0167     [P,~] = polyfit(x,y,1);
0168 <span class="comment">%     figure(1); hold on;</span>
0169 <span class="comment">%     plot(x,polyval(P,x),'g-')</span>
0170     V.m = P(1);
0171     V.b = P(2);
0172 <span class="keyword">else</span>
0173     [P,~] = polyfit(y,x,1);
0174 <span class="comment">%     figure(1); hold on;</span>
0175 <span class="comment">%     plot(polyval(P,y),y,'g-')</span>
0176     V.m = 1/P(1);           <span class="comment">%Reformat slope and intercept in terms of y= fn(x) rather than x = fn(y)</span>
0177     V.b = -P(2)/P(1);
0178 <span class="keyword">end</span>
0179 
0180 <span class="comment">%Former method commented out</span>
0181 <span class="comment">% whichstats = {'tstat','yhat'};</span>
0182 <span class="comment">% stats = regstats(y,x,'linear', whichstats);</span>
0183 <span class="comment">%</span>
0184 <span class="comment">% % mean cross-section line slope and intercept</span>
0185 <span class="comment">% V.m = stats.tstat.beta(2);</span>
0186 <span class="comment">% V.b = stats.tstat.beta(1);</span>
0187 
0188 clear x y stats whichstats zi
0189 
0190 <span class="comment">%% Map ensembles to mean c-s line</span>
0191 <span class="comment">% Determine the point (mapped ensemble point) where the equation of the</span>
0192 <span class="comment">% mean cross-section line intercepts a line perpendicular to the mean</span>
0193 <span class="comment">% cross-section line passing through an ensemble from an individual</span>
0194 <span class="comment">% transect (see notes for equation derivation)</span>
0195 
0196 <span class="keyword">for</span> zi = 1:z
0197     A(zi).Comp.xm = ((A(zi).Comp.xUTM-V.m.*V.b+V.m.*A(zi).Comp.yUTM)<span class="keyword">...</span>
0198         ./(V.m.^2+1));
0199     A(zi).Comp.ym = ((V.b+V.m.*A(zi).Comp.xUTM+V.m.^2.*A(zi).Comp.yUTM)<span class="keyword">...</span>
0200         ./(V.m.^2+1));
0201 <span class="keyword">end</span>
0202 
0203 <span class="comment">%Plot data to check</span>
0204 xensall = [];
0205 yensall = [];
0206 <span class="keyword">for</span> zi = 1:z
0207 <span class="comment">%   plot(A(zi).Comp.xm,A(zi).Comp.ym,'b.')</span>
0208   xensall = [xensall; A(zi).Comp.xm];
0209   yensall = [yensall; A(zi).Comp.ym];
0210 <span class="keyword">end</span>
0211 <span class="comment">% % plot(A(3).Comp.xm,A(3).Comp.ym,'xg')</span>
0212 <span class="comment">% % plot(A(4).Comp.xm,A(4).Comp.ym,'oy')</span>
0213 <span class="comment">% xlabel('UTM Easting (m)')</span>
0214 <span class="comment">% ylabel('UTM Northing (m)')</span>
0215 <span class="comment">% box on</span>
0216 <span class="comment">% grid on</span>
0217 <span class="comment">% %Plot a legend in Figure 1</span>
0218 <span class="comment">% %figure(1); hold on</span>
0219 <span class="comment">% %legend('Shoreline','GPS(corr)','GPS(raw)','Best Fit','Trans 1</span>
0220 <span class="comment">% %(mapped)','Other Trans (mapped)')</span>
0221 
0222 <span class="comment">%Compute the median distance between mapped points</span>
0223 Dmat = [xensall yensall];
0224 <span class="keyword">if</span> xrng &gt; yrng
0225     Dmat = sortrows(Dmat,1);
0226 <span class="keyword">else</span>
0227     Dmat = sortrows(Dmat,2);
0228 <span class="keyword">end</span>
0229 dxall = diff(Dmat(:,1));
0230 dyall = diff(Dmat(:,2));
0231 densall = sqrt(dxall.^2 + dyall.^2);
0232 V.meddens = median(densall);
0233 V.stddens = std(densall);
0234 disp([<span class="string">'Median Spacing Between Mapped Ensembles = '</span> num2str(V.meddens) <span class="string">' m'</span>])
0235 disp([<span class="string">'Standard Deviation of Spacing Between Mapped Ensembles = '</span> num2str(V.stddens) <span class="string">' m'</span>])
0236 disp([<span class="string">'Recommended Grid Node Spacing &gt; '</span> num2str(V.meddens + V.stddens) <span class="string">' m'</span>])
0237 
0238 <span class="comment">%Display in message box for compiled version</span>
0239 msg_string = {[<span class="string">'Median Spacing Between Mapped Ensembles = '</span> num2str(V.meddens) <span class="string">' m'</span>],<span class="keyword">...</span>
0240     [<span class="string">'Standard Deviation of Spacing Between Mapped Ensembles = '</span> num2str(V.stddens) <span class="string">' m'</span>],<span class="keyword">...</span>
0241     [<span class="string">'Recommended Grid Node Spacing &gt; '</span> num2str(V.meddens + V.stddens) <span class="string">' m'</span>]};
0242 msgbox(msg_string,<span class="string">'VMT Grid Node Spacing'</span>,<span class="string">'help'</span>,<span class="string">'replace'</span>);
0243 
0244 <span class="comment">%%Save the shorepath</span>
0245 <span class="keyword">if</span> setends
0246     <span class="keyword">if</span> exist(<span class="string">'LastDir.mat'</span>,<span class="string">'dir'</span>)
0247         save(<span class="string">'LastDir.mat'</span>,<span class="string">'endspath'</span>,<span class="string">'-append'</span>)
0248     <span class="keyword">else</span>
0249         save(<span class="string">'LastDir.mat'</span>,<span class="string">'endspath'</span>)
0250     <span class="keyword">end</span>
0251 <span class="keyword">end</span>
0252 
0253 <span class="comment">%==========================================================================</span>
0254 <a name="_sub2" href="#_subfunctions" class="code">function [A,V] = VMT_GridData2MeanXS(z,A,V)</a>
0255 
0256 <span class="comment">%This routine generates a uniformly spaced grid for the mean cross section and</span>
0257 <span class="comment">%maps (interpolates) individual transects to this grid.</span>
0258 
0259 <span class="comment">%(adapted from code by J. Czuba)</span>
0260 
0261 <span class="comment">%P.R. Jackson, USGS, 12-9-08</span>
0262 
0263 <span class="comment">%% User Input</span>
0264 
0265 xgdspc = A(1).hgns; <span class="comment">%Horizontal Grid node spacing in meters  (vertical grid spacing is set by bins)</span>
0266 <span class="comment">% if 0</span>
0267 <span class="comment">%     xgdspc = V.meddens + V.stddens;  %Auto method should include 67% of the values</span>
0268 <span class="comment">%     disp(['X Grid Node Auto Spacing = ' num2str(xgdspc) ' m'])</span>
0269 <span class="comment">% end</span>
0270 
0271 <span class="comment">%% Determine uniform mean c-s grid for vector interpolating</span>
0272 <span class="comment">% Determine the end points of the mean cross-section line</span>
0273 <span class="comment">% Initialize variable with mid range value</span>
0274 V.xe = mean(A(1).Comp.xm);
0275 V.ys = mean(A(1).Comp.ym);
0276 V.xw = mean(A(1).Comp.xm);
0277 V.yn = mean(A(1).Comp.ym);
0278 
0279 <span class="keyword">for</span> zi = 1:z
0280     V.xe = max(max(A(zi).Comp.xm),V.xe);
0281     V.ys = min(min(A(zi).Comp.ym),V.ys);
0282     V.xw = min(min(A(zi).Comp.xm),V.xw);
0283     V.yn = max(max(A(zi).Comp.ym),V.yn);
0284 <span class="keyword">end</span>
0285 
0286 <span class="comment">% Determine the distance between the mean cross-section endpoints</span>
0287 V.dx = V.xe-V.xw;
0288 V.dy = V.yn-V.ys;
0289 
0290 V.dl = sqrt(V.dx^2+V.dy^2);
0291 
0292 <span class="comment">% Determine mean cross-section velocity vector grid</span>
0293 V.mcsDist = linspace(0,V.dl,floor(V.dl/xgdspc));                                  <span class="comment">%%linspace(0,V.dl,V.dl); Changed to make it user selectable (PRJ, 12-12-08)</span>
0294 V.mcsDepth = A(1).Wat.binDepth(:,1);
0295 [V.mcsDist V.mcsDepth] = meshgrid(V.mcsDist,V.mcsDepth);
0296 
0297 <span class="comment">% Determine the angle the mean cross-section makes with the</span>
0298 <span class="comment">% x-axis (longitude)</span>
0299 <span class="comment">% Plot mean cross-section line</span>
0300 <span class="keyword">if</span> V.m &gt;= 0
0301     V.theta = atand(V.dy./V.dx);
0302     
0303 <span class="comment">%     figure(1); hold on</span>
0304 <span class="comment">%     plot([V.xe V.xw],[V.yn V.ys],'ks'); hold on</span>
0305     
0306     V.mcsX = V.xe - V.mcsDist(1,:).*cosd(V.theta);            <span class="comment">%</span>
0307     V.mcsY = V.yn - V.mcsDist(1,:).*sind(V.theta); 
0308     
0309 <span class="comment">%     if V.mfd &gt;= 270 | V.mfd &lt; 90 %Flow to the north  %This code was an attempt to auto detect left bank--did'nt work so removed.</span>
0310 <span class="comment">%         V.mcsX = V.xw+V.mcsDist(1,:).*cosd(V.theta);            %</span>
0311 <span class="comment">%         V.mcsY = V.ys+V.mcsDist(1,:).*sind(V.theta);</span>
0312 <span class="comment">%</span>
0313 <span class="comment">%     elseif V.mfd &gt;= 90 &amp; V.mfd &lt; 270 %Flow to the south</span>
0314 <span class="comment">%         V.mcsX = V.xe-V.mcsDist(1,:).*cosd(V.theta);            %</span>
0315 <span class="comment">%         V.mcsY = V.yn-V.mcsDist(1,:).*sind(V.theta);</span>
0316 <span class="comment">%     end%</span>
0317     
0318 <span class="comment">%     plot(V.mcsX,V.mcsY,'k+'); hold on</span>
0319 <span class="comment">%     plot(V.mcsX(1),V.mcsY(1),'y*'); hold on</span>
0320 
0321 <span class="keyword">elseif</span> V.m &lt; 0
0322     V.theta = -atand(V.dy./V.dx);
0323     <span class="comment">%V.theta = atand(V.dy./V.dx); %Changed 9-28-10, PRJ (theta needs to be</span>
0324     <span class="comment">%negative--changed back to original)</span>
0325     
0326 <span class="comment">%     figure(1); hold on</span>
0327 <span class="comment">%     plot([V.xe V.xw],[V.ys V.yn],'ks'); hold on</span>
0328     
0329     V.mcsX = V.xe - V.mcsDist(1,:).*cosd(V.theta);
0330     V.mcsY = V.ys - V.mcsDist(1,:).*sind(V.theta);
0331     <span class="comment">%V.mcsY = V.ys+V.mcsDist(1,:).*sind(V.theta);  %Changed 9-28-10, PRJ</span>
0332     
0333 <span class="comment">%     if V.mfd &gt;= 270 | V.mfd &lt; 90 %Flow to the north</span>
0334 <span class="comment">%         V.mcsX = V.xw+V.mcsDist(1,:).*cosd(V.theta);            %</span>
0335 <span class="comment">%         V.mcsY = V.yn+V.mcsDist(1,:).*sind(V.theta);</span>
0336 <span class="comment">%</span>
0337 <span class="comment">%     elseif V.mfd &gt;= 90 &amp; V.mfd &lt; 270 %Flow to the south</span>
0338 <span class="comment">%         V.mcsX = V.xe-V.mcsDist(1,:).*cosd(V.theta);</span>
0339 <span class="comment">%         V.mcsY = V.ys-V.mcsDist(1,:).*sind(V.theta);</span>
0340 <span class="comment">%     end%</span>
0341    
0342 <span class="comment">%     plot(V.mcsX,V.mcsY,'k+'); hold on</span>
0343 <span class="comment">%     plot(V.mcsX(1),V.mcsY(1),'y*'); hold on</span>
0344 <span class="keyword">end</span>
0345 
0346 V.mcsX = meshgrid(V.mcsX,V.mcsDepth(:,1));
0347 V.mcsY = meshgrid(V.mcsY,V.mcsDepth(:,1));
0348 <span class="comment">% figure(1); set(gca,'DataAspectRatio',[1 1 1],'PlotBoxAspectRatio',[1 1 1])</span>
0349 clear zi
0350 
0351 <span class="comment">% Format the ticks for UTM and allow zooming and panning</span>
0352 <span class="comment">% figure(1);</span>
0353 <span class="comment">% ticks_format('%6.0f','%8.0f'); %formats the ticks for UTM</span>
0354 <span class="comment">% hdlzm_fig1 = zoom;</span>
0355 <span class="comment">% set(hdlzm_fig1,'ActionPostCallback',@mypostcallback_zoom);</span>
0356 <span class="comment">% set(hdlzm_fig1,'Enable','on');</span>
0357 <span class="comment">% hdlpn_fig1 = pan;</span>
0358 <span class="comment">% set(hdlpn_fig1,'ActionPostCallback',@mypostcallback_pan);</span>
0359 <span class="comment">% set(hdlpn_fig1,'Enable','on');</span>
0360 
0361 
0362 <span class="comment">%% Determine location of mapped ensemble points for interpolating</span>
0363 <span class="comment">% For all transects</span>
0364 
0365 <span class="comment">%A = VMT_DxDyfromLB(z,A,V); %Computes dx and dy</span>
0366 
0367 <span class="keyword">for</span> zi = 1:z
0368     <span class="comment">% Determine if the mean cross-section line trends NW-SE or SW-NE</span>
0369     <span class="comment">% Determine the distance in radians from the left bank mean</span>
0370     <span class="comment">% cross-section point to the mapped ensemble point for an individual</span>
0371     <span class="comment">% transect</span>
0372     A(zi).Comp.dx = abs(V.xe-A(zi).Comp.xm);  <span class="comment">%This assumes the easternmost bank is the left bank--changed PRJ 1-21-09 (now uses VMT_DxDyfromLB--not working 2/1/09)</span>
0373 
0374     <span class="keyword">if</span> V.m &gt; 0
0375         A(zi).Comp.dy = abs(V.yn-A(zi).Comp.ym);
0376     <span class="keyword">elseif</span> V.m &lt; 0
0377         A(zi).Comp.dy = abs(V.ys-A(zi).Comp.ym);
0378     <span class="keyword">end</span>
0379 
0380     <span class="comment">% Determine the distance in meters from the left bank mean</span>
0381     <span class="comment">% cross-section point to the mapped ensemble point for an individual</span>
0382     <span class="comment">% transect</span>
0383     A(zi).Comp.dl = sqrt(A(zi).Comp.dx.^2+A(zi).Comp.dy.^2);
0384 
0385 <span class="comment">%     % Sort vectors by dl</span>
0386 <span class="comment">%     A(zi).Comp.dlsort = sort(A(zi).Comp.dl,'ascend');</span>
0387 <span class="comment">%</span>
0388 <span class="comment">%     % Map indices</span>
0389 <span class="comment">%     for i = 1:A(zi).Sup.noe</span>
0390 <span class="comment">%         for k = 1 : A(zi).Sup.noe</span>
0391 <span class="comment">%             if A(zi).Comp.dlsort(i,1) == A(zi).Comp.dl(k,1)</span>
0392 <span class="comment">%                 A(zi).Comp.vecmap(i,1) = k;</span>
0393 <span class="comment">%             end</span>
0394 <span class="comment">%         end</span>
0395 <span class="comment">%     end</span>
0396     
0397     <span class="comment">% Sort vectors by dl</span>
0398     [A(zi).Comp.dlsort,A(zi).Comp.vecmap] = sort(A(zi).Comp.dl,<span class="string">'ascend'</span>);
0399 
0400     <span class="comment">% GPS position fix</span>
0401     <span class="comment">% if distances from the left bank are the same for two ensembles the</span>
0402     <span class="comment">% the position of the right most ensemble is interpolated from adjacent</span>
0403     <span class="comment">% ensembles</span>
0404     <span class="comment">% check for repeat values of distance</span>
0405     sbt(:,1) = diff(A(zi).Comp.dlsort);
0406     chk(1,1) = 1;
0407     chk(2:A(zi).Sup.noe,1) = sbt(1:<span class="keyword">end</span>,1);
0408 
0409     <span class="comment">% identify repeat values</span>
0410     A(zi).Comp.sd = (chk==0) &gt; 0;
0411 
0412     <span class="comment">% if repeat values exist interpolate distances from adjacent ensembles</span>
0413     <span class="keyword">if</span> sum(A(zi).Comp.sd) &gt; 0
0414 
0415         <span class="comment">% bracket repeat sections</span>
0416         [I,~] = ind2sub(size(A(zi).Comp.sd),find(A(zi).Comp.sd==1));
0417         df = diff(I);
0418         nbrk = sum(df&gt;1)+1;
0419         [I2,~] = ind2sub(size(df),find(df&gt;1));
0420 
0421         bg(1) = I(1)-1;
0422         <span class="keyword">for</span> n = 2 : nbrk
0423             bg(n) = I(I2(n-1)+1)-1;
0424         <span class="keyword">end</span>
0425 
0426         <span class="keyword">for</span> n = 1:nbrk -1
0427             ed(n) = I(I2(n))+1;
0428         <span class="keyword">end</span>
0429         ed(nbrk) = I(end)+1;
0430 
0431         <span class="comment">% interpolate repeat values</span>
0432         A(zi).Comp.dlsortgpsfix = A(zi).Comp.dlsort;
0433 
0434         <span class="keyword">for</span> i = 1 : nbrk
0435             <span class="keyword">for</span> j = bg(i)+1 : ed(i)-1
0436                 <span class="comment">% interpolate</span>
0437                 <span class="keyword">if</span> bg(i) &gt; 0 &amp;&amp; ed(i) &lt; length(A(zi).Nav.lat_deg)
0438 
0439                     den=(ed(i)-bg(i));
0440                     num2=j-bg(i);
0441                     num1=ed(i)-j;
0442 
0443                     A(zi).Comp.dlsortgpsfix(j,1)=<span class="keyword">...</span>
0444                         (num1/den).*A(zi).Comp.dlsort(bg(i))+<span class="keyword">...</span>
0445                         (num2/den).*A(zi).Comp.dlsort(ed(i));
0446 
0447                 <span class="keyword">end</span>
0448                 
0449                 <span class="comment">% extrapolate end</span>
0450                 <span class="keyword">if</span> ed(i) &gt; length(A(zi).Nav.lat_deg)
0451                    
0452                     numex = ed(i) - length(A(zi).Nav.lat_deg);
0453                     
0454                     A(zi).Comp.dlsortgpsfix(j,1)=numex.*<span class="keyword">...</span>
0455                         (A(zi).Comp.dlsort(bg(i))-<span class="keyword">...</span>
0456                         A(zi).Comp.dlsort(bg(i)-1))+<span class="keyword">...</span>
0457                         A(zi).Comp.dlsort(bg(i));
0458                     
0459                 <span class="keyword">end</span>               
0460             <span class="keyword">end</span>
0461         <span class="keyword">end</span>
0462     <span class="keyword">else</span>
0463         A(zi).Comp.dlsortgpsfix = A(zi).Comp.dlsort;
0464     <span class="keyword">end</span> <span class="comment">% IF sum(A(zi).Comp.sd) &gt; 0</span>
0465     
0466     <span class="comment">% Determine velocity vector grid for individual transects</span>
0467     [A(zi).Comp.itDist A(zi).Comp.itDepth] = <span class="keyword">...</span>
0468         meshgrid(A(zi).Comp.dlsortgpsfix,A(zi).Wat.binDepth(:,1));
0469     
0470     clear I I2 J J2 bg chk df ed i j nbrk sbt xUTM yUTM n zi<span class="keyword">...</span>
0471         den num2 num1 numex
0472     
0473 <span class="keyword">end</span>
0474 
0475 clear zi i k check
0476 
0477 <span class="comment">%% Interpolate individual transects onto uniform mean c-s grid</span>
0478 <span class="comment">% Fill in uniform grid based on individual transects mapped onto the mean</span>
0479 <span class="comment">% cross-section by interpolating between adjacent points</span>
0480 
0481 <span class="comment">%ZI = interp2(X,Y,Z,XI,YI)</span>
0482 <span class="keyword">for</span> zi = 1:z
0483     A(zi).Comp.mcsBack = <span class="keyword">...</span>
0484         interp2(A(zi).Comp.itDist, <span class="keyword">...</span>
0485                 A(zi).Comp.itDepth, <span class="keyword">...</span>
0486                 A(zi).Clean.bs(:,A(zi).Comp.vecmap), <span class="keyword">...</span>
0487                 V.mcsDist, <span class="keyword">...</span>
0488                 V.mcsDepth);
0489     A(zi).Comp.mcsBack(A(zi).Comp.mcsBack&gt;=255) = NaN;
0490     
0491     <span class="comment">%A(zi).Comp.mcsDir = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...</span>
0492         <span class="comment">%A(zi).Clean.vDir(:,A(zi).Comp.vecmap), V.mcsDist, V.mcsDepth); %This interpolation scheme has issues when interpolating in a flow due north (0,360 interpolate to 180)</span>
0493     
0494     <span class="comment">% For direction, must convert degrees to radians, take the sin of the</span>
0495     <span class="comment">% radians, and then interpolate.  Following interpolation, convert</span>
0496     <span class="comment">% radians back to degrees. (PRJ, 9-28-10)  ALSO BAD NEAR 180</span>
0497     <span class="comment">%A(zi).Comp.mcsDir = 180/pi*(interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...</span>
0498         <span class="comment">%sin(pi/180*(A(zi).Clean.vDir(:,A(zi).Comp.vecmap))), V.mcsDist, V.mcsDepth));</span>
0499     
0500 <span class="comment">%     A(zi).Comp.mcsMag = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...</span>
0501 <span class="comment">%         A(zi).Clean.vMag(:,A(zi).Comp.vecmap), V.mcsDist, V.mcsDepth);</span>
0502 <span class="comment">%         (Recomputed from north and east components (PRJ, 3-21-11)</span>
0503     
0504     
0505     A(zi).Comp.mcsEast = <span class="keyword">...</span>
0506         interp2(A(zi).Comp.itDist, <span class="keyword">...</span>
0507                 A(zi).Comp.itDepth, <span class="keyword">...</span>
0508                 A(zi).Clean.vEast(:,A(zi).Comp.vecmap), <span class="keyword">...</span>
0509                 V.mcsDist, <span class="keyword">...</span>
0510                 V.mcsDepth);
0511     A(zi).Comp.mcsNorth = <span class="keyword">...</span>
0512         interp2(A(zi).Comp.itDist, <span class="keyword">...</span>
0513                 A(zi).Comp.itDepth, <span class="keyword">...</span>
0514                 A(zi).Clean.vNorth(:,A(zi).Comp.vecmap), <span class="keyword">...</span>
0515                 V.mcsDist, <span class="keyword">...</span>
0516                 V.mcsDepth);
0517     A(zi).Comp.mcsVert = <span class="keyword">...</span>
0518         interp2(A(zi).Comp.itDist, <span class="keyword">...</span>
0519                 A(zi).Comp.itDepth, <span class="keyword">...</span>
0520                 A(zi).Clean.vVert(:,A(zi).Comp.vecmap), <span class="keyword">...</span>
0521                 V.mcsDist, <span class="keyword">...</span>
0522                 V.mcsDepth);
0523     
0524     <span class="comment">%Compute magnitude</span>
0525     A(zi).Comp.mcsMag = sqrt(A(zi).Comp.mcsEast.^2 + A(zi).Comp.mcsNorth.^2);
0526     
0527     <span class="comment">%For direction, compute from the velocity components</span>
0528     A(zi).Comp.mcsDir = 90 - (atan2(A(zi).Comp.mcsNorth, A(zi).Comp.mcsEast))*180/pi; <span class="comment">%Compute the atan from the velocity componentes, convert to radians, and rotate to north axis</span>
0529     qindx = find(A(zi).Comp.mcsDir &lt; 0);
0530     <span class="keyword">if</span> ~isempty(qindx)
0531         A(zi).Comp.mcsDir(qindx) = A(zi).Comp.mcsDir(qindx) + 360;  <span class="comment">%Must add 360 deg to Quadrant 4 values as they are negative angles from the +y axis</span>
0532     <span class="keyword">end</span>
0533     
0534     A(zi).Comp.mcsBed = <span class="keyword">...</span>
0535         interp1(A(zi).Comp.itDist(1,:), <span class="keyword">...</span>
0536                 nanmean(A(zi).Nav.depth(A(zi).Comp.vecmap,:),2), <span class="keyword">...</span>
0537                 V.mcsDist(1,:));
0538 <span class="keyword">end</span>
0539 
0540 clear zi
0541 
0542 <span class="comment">%% Embedded functions</span>
0543 <a name="_sub3" href="#_subfunctions" class="code">function mypostcallback_zoom(obj,evd)</a>
0544 ticks_format(<span class="string">'%6.0f'</span>,<span class="string">'%8.0f'</span>); <span class="comment">%formats the ticks for UTM (when zooming)</span>
0545 
0546 <a name="_sub4" href="#_subfunctions" class="code">function mypostcallback_pan(obj,evd)</a>
0547 ticks_format(<span class="string">'%6.0f'</span>,<span class="string">'%8.0f'</span>); <span class="comment">%formats the ticks for UTM (when panning)</span>
0548 
0549 <span class="comment">%==========================================================================</span>
0550 <a name="_sub5" href="#_subfunctions" class="code">function [A,V] = VMT_CompMeanXS_old(z,A,V)</a>
0551 <span class="comment">% Compute the mean cross section data from individual transects that have</span>
0552 <span class="comment">% been previously mapped to a common grid.</span>
0553 
0554 <span class="comment">%(adapted from code by J. Czuba)</span>
0555 
0556 <span class="comment">%P.R. Jackson, USGS, 12-9-08</span>
0557 
0558 <span class="comment">%% Average mapped mean cross-sections from individual transects together</span>
0559 <span class="comment">% Assign mapped uniform grid vectors to the same array for averaging</span>
0560 Back  = zeros([size(A(1).Comp.mcsBack) z]);
0561 Dir   = zeros([size(A(1).Comp.mcsDir) z]);
0562 Mag   = zeros([size(A(1).Comp.mcsMag) z]);
0563 East  = zeros([size(A(1).Comp.mcsEast) z]);
0564 North = zeros([size(A(1).Comp.mcsNorth) z]);
0565 Vert  = zeros([size(A(1).Comp.mcsVert) z]);
0566 Bed   = zeros([size(A(1).Comp.mcsBed) z]);
0567 <span class="keyword">for</span> zi = 1 : z
0568     Back(:,:,zi)  = A(zi).Comp.mcsBack;
0569     Dir(:,:,zi)   = A(zi).Comp.mcsDir;
0570     Mag(:,:,zi)   = A(zi).Comp.mcsMag;
0571     East(:,:,zi)  = A(zi).Comp.mcsEast;
0572     North(:,:,zi) = A(zi).Comp.mcsNorth;
0573     Vert(:,:,zi)  = A(zi).Comp.mcsVert;
0574     Bed(:,:,zi)   = A(zi).Comp.mcsBed;
0575 <span class="keyword">end</span>
0576 
0577 <span class="comment">% numavg = nansum(~isnan(Mag),3);</span>
0578 <span class="comment">% numavg(numavg==0) = NaN;</span>
0579 <span class="comment">% enscnt = nanmean(numavg,1);</span>
0580 <span class="comment">% [I,J] = ind2sub(size(enscnt),find(enscnt&gt;=1));  %Changed to &gt;= 1 PRJ 12-10-08  (uses data even if only one measurement)</span>
0581 
0582 Backone = Back;
0583 Magone  = Mag;
0584 Vertone = Vert;
0585 Bedone  = Bed;
0586 
0587 Backone(~isnan(Back)) = 1;
0588 Magone(~isnan(Mag))   = 1;
0589 Vertone(~isnan(Vert)) = 1;
0590 Bedone(~isnan(Bed))   = 1;
0591 
0592 V.countBack = nansum(Backone,3);
0593 V.countMag  = nansum(Magone,3);
0594 V.countVert = nansum(Vertone,3);
0595 V.countBed  = nansum(Bedone,3);
0596 
0597 V.countBack(V.countBack==0) = NaN;
0598 V.countMag(V.countMag==0)   = NaN;
0599 V.countVert(V.countVert==0) = NaN;
0600 V.countBed(V.countBed==0)   = NaN;
0601 
0602 <span class="comment">% Average mapped mean cross-sections from individual transects together</span>
0603 V.mcsBack  = nanmean(Back,3);
0604 <span class="comment">% V.mcsDir   = nanmean(Dir,3);  % Will not average correctly in all cases due to 0-360</span>
0605 <span class="comment">%wrapping (PRJ, 9-29-10)</span>
0606 <span class="comment">% V.mcsMag   = nanmean(Mag,3);  %Mag recomputed from north, east, up components(PRJ, 3-21-11)</span>
0607 V.mcsEast  = nanmean(East,3);
0608 V.mcsNorth = nanmean(North,3);
0609 V.mcsVert  = nanmean(Vert,3);
0610 
0611 <span class="comment">% Average Magnitude</span>
0612 V.mcsMag = sqrt(V.mcsEast.^2 + V.mcsNorth.^2 + V.mcsVert.^2);
0613 
0614 <span class="comment">% Average the flow direction</span>
0615 V.mcsDir = 90 - (atan2(V.mcsNorth, V.mcsEast))*180/pi; <span class="comment">%Compute the atan from the velocity componentes, convert to radians, and rotate to north axis</span>
0616 qindx = find(V.mcsDir &lt; 0);
0617 <span class="keyword">if</span> ~isempty(qindx)
0618     V.mcsDir(qindx) = V.mcsDir(qindx) + 360;  <span class="comment">%Must add 360 deg to Quadrant 4 values as they are negative angles from the +y axis</span>
0619 <span class="keyword">end</span>
0620 
0621 V.mcsBed = nanmean(Bed,3);
0622 
0623 <span class="comment">% Compute the Bed Elevation in meters</span>
0624 disp([<span class="string">'Assigned Water Surface Elevation (WSE; in meters) = '</span> num2str(A(1).wse)])
0625 V.mcsBedElev = A(1).wse - V.mcsBed;
0626 
0627 <span class="comment">%==========================================================================</span>
0628 <a name="_sub6" href="#_subfunctions" class="code">function [A,V] = VMT_CompMeanXS_UVW(z,A,V)</a>
0629 
0630 <span class="comment">%This routine computes the mean cross section velocity components (UVW)</span>
0631 <span class="comment">%from individual transects that have been previously mapped to a common grid and averaged.</span>
0632 
0633 <span class="comment">%(adapted from code by J. Czuba)</span>
0634 
0635 <span class="comment">%P.R. Jackson, USGS, 12-9-08</span>
0636 
0637 
0638 <span class="comment">%% Rotate velocities into u, v, and w components</span>
0639 <span class="comment">% Determine the direction of streamwise velocity (u)</span>
0640 V.phi = 180 - V.theta;  <span class="comment">%Taken as perpendicular to the mean XS</span>
0641 
0642 <span class="comment">% Determine the deviation of a vector from streamwise velocity</span>
0643 V.psi = V.phi - V.mcsDir;
0644 
0645 <span class="comment">% Determine streamwise (u), transverse (v), and vertical (w) velocities</span>
0646 V.u = cosd(V.psi).*V.mcsMag;
0647 V.v = sind(V.psi).*V.mcsMag;
0648 V.w = V.mcsVert;
0649 
0650 <span class="keyword">for</span> zi = 1:z
0651     A(zi).Comp.u = cosd(V.psi).*A(zi).Comp.mcsMag;
0652     A(zi).Comp.v = sind(V.psi).*A(zi).Comp.mcsMag;
0653     A(zi).Comp.w = A(zi).Comp.mcsVert;
0654 <span class="keyword">end</span>
0655 
0656 <span class="comment">%==========================================================================</span>
0657 <a name="_sub7" href="#_subfunctions" class="code">function [A,V] = VMT_CompMeanXS_PriSec(z,A,V)</a>
0658 
0659 <span class="comment">%This routine computes the mean cross section velocity components (Primary and secondary)</span>
0660 <span class="comment">%from individual transects that have been previously mapped to a common grid and averaged.</span>
0661 <span class="comment">%The Primary velocity is defined as the component of the flow in the direction of the discharge</span>
0662 <span class="comment">%(i.e. rotated from the streamwise direction so the secrondary discharge is</span>
0663 <span class="comment">%zero).</span>
0664 
0665 <span class="comment">% This is referred to as the &quot;zero net cross-stream discharge definition&quot;</span>
0666 <span class="comment">% (see Lane et al. 2000, Hydrological Processes 14, 2047-2071)</span>
0667 
0668 <span class="comment">%(adapted from code by J. Czuba)</span>
0669 
0670 <span class="comment">%P.R. Jackson, USGS, 12-9-08</span>
0671 
0672 <span class="comment">%% Rotate velocities into p and s components for the mean transect</span>
0673 <span class="comment">% calculate dy and dz for each meaurement point</span>
0674 dy = mean(diff(V.mcsDist(1,:)));  <span class="comment">% m</span>
0675 dz = mean(diff(V.mcsDepth(:,1))); <span class="comment">% m</span>
0676 dydz = dy.*dz;
0677 
0678 <span class="comment">% calculate the bit of discharge for each imaginary cell around the</span>
0679 <span class="comment">% velocity point</span>
0680 qyi = V.v.*dydz; <span class="comment">% cm*m^2/s</span>
0681 qxi = V.u.*dydz; <span class="comment">% cm*m^2/s</span>
0682 <span class="comment">% qyi = V.v.*dy.*dz;%cm*m^2/s</span>
0683 <span class="comment">% qxi = V.u.*dy.*dz;%cm*m^2/s</span>
0684 
0685 <span class="comment">% sum the streamwise and transverse Q and calculate the angle of the</span>
0686 <span class="comment">% cross section</span>
0687 V.Qy = nansum(nansum(qyi));<span class="comment">%cm*m^2/s</span>
0688 V.Qx = nansum(nansum(qxi));<span class="comment">%cm*m^2/s</span>
0689 <span class="comment">% V.Qy = nansum(qyi(:)); % cm*m^2/s % PDF?</span>
0690 <span class="comment">% V.Qx = nansum(qxi(:)); % cm*m^2/s</span>
0691 
0692 V.alphasp = atand(V.Qy./V.Qx);
0693 V.phisp   = V.phi - V.alphasp;
0694 
0695 <span class="comment">% Rotate the velocities so that Qy is effectively zero</span>
0696 qpi =  qxi.*cosd(V.alphasp) + qyi.*sind(V.alphasp);
0697 qsi = -qxi.*sind(V.alphasp) + qyi.*cosd(V.alphasp);
0698 <span class="comment">% R =  [ cosd(V.alphasp) sind(V.alphasp)   % PDF? Depends on dimensions</span>
0699 <span class="comment">%       -sind(V.alphasp) cosd(V.alphasp)];</span>
0700 
0701 V.Qp = nansum(nansum(qpi));<span class="comment">%cm*m^2/s</span>
0702 V.Qs = nansum(nansum(qsi));<span class="comment">%cm*m^2/s</span>
0703 <span class="comment">% V.Qp = nansum(qpi(:)));%cm*m^2/s % PDF?</span>
0704 <span class="comment">% V.Qs = nansum(qsi(:)));%cm*m^2/s</span>
0705 disp([<span class="string">'Secondary Discharge after Rotation (ZSD definition; m^3/s) = '</span> num2str(V.Qs/100)])
0706 
0707 V.vp = qpi./dydz; <span class="comment">% cm/s</span>
0708 V.vs = qsi./dydz; <span class="comment">% cm/s</span>
0709 <span class="comment">% V.vp = qpi./(dy.*dz);%cm/s</span>
0710 <span class="comment">% V.vs = qsi./(dy.*dz);%cm/s</span>
0711 
0712 <span class="comment">%% Transform each individual transect</span>
0713 
0714 <span class="keyword">for</span> zi = 1 : z
0715     <span class="comment">% calculate the bit of discharge for each imaginary cell around the</span>
0716     <span class="comment">% velocity point</span>
0717     A(zi).Comp.qyi = A(zi).Comp.v.*dydz; <span class="comment">% cm*m^2/s</span>
0718     A(zi).Comp.qxi = A(zi).Comp.u.*dydz; <span class="comment">% cm*m^2/s</span>
0719 <span class="comment">%     A(zi).Comp.qyi=A(zi).Comp.v.*dy.*dz;%cm*m^2/s</span>
0720 <span class="comment">%     A(zi).Comp.qxi=A(zi).Comp.u.*dy.*dz;%cm*m^2/s</span>
0721     
0722     <span class="comment">% rotate the velocities so that Qy is effcetively zero</span>
0723     A(zi).Comp.qpi =  A(zi).Comp.qxi.*cosd(V.alphasp) + A(zi).Comp.qyi.*sind(V.alphasp);
0724     A(zi).Comp.qsi = -A(zi).Comp.qxi.*sind(V.alphasp) + A(zi).Comp.qyi.*cosd(V.alphasp);
0725     
0726     A(zi).Comp.Qp=nansum(nansum(A(zi).Comp.qpi));<span class="comment">%cm*m^2/s</span>
0727     A(zi).Comp.Qs=nansum(nansum(A(zi).Comp.qsi));<span class="comment">%cm*m^2/s</span>
0728     <span class="comment">% A(zi).Comp.Qp = nansum(A(zi).Comp.qpi(:)); % cm*m^2/s  % PDF?</span>
0729     <span class="comment">% A(zi).Comp.Qs = nansum(A(zi).Comp.qsi(:)); % cm*m^2/s</span>
0730     
0731     A(zi).Comp.vp = A(zi).Comp.qpi./dydz; <span class="comment">% cm/s</span>
0732     A(zi).Comp.vs  =A(zi).Comp.qsi./dydz; <span class="comment">% cm/s</span>
0733 <span class="comment">%     A(zi).Comp.vp=A(zi).Comp.qpi./(dy.*dz);%cm/s</span>
0734 <span class="comment">%     A(zi).Comp.vs=A(zi).Comp.qsi./(dy.*dz);%cm/s</span>
0735 <span class="keyword">end</span>
0736 
0737 
0738 <span class="comment">%% Determine velocity deviations from the p direction</span>
0739 
0740 V.mcsDirDevp = V.phisp - V.mcsDir;
0741 
0742 <span class="keyword">for</span> zi = 1:z
0743     A(zi).Comp.mcsDirDevp = V.phisp - A(zi).Comp.mcsDir;
0744 <span class="keyword">end</span>
0745 
0746 <span class="comment">%==========================================================================</span>
0747 <a name="_sub8" href="#_subfunctions" class="code">function [V] = VMT_RozovskiiV2(V,A)</a>
0748 
0749 <span class="comment">% Computes the frame of reference transpositon as described in Kenworthy</span>
0750 <span class="comment">% and Rhoads (1998) ESPL using a Rozovskii type analysis.</span>
0751 
0752 <span class="comment">% Notes:</span>
0753 <span class="comment">%     -The depth averaging currently linearly interpolates to the bed,</span>
0754 <span class="comment">%      however we may want some other approach such as log law, etc.</span>
0755 <span class="comment">%     -I extrapolate the velocity at the near surface bin to the water</span>
0756 <span class="comment">%      surface for the depth averaging (ie, BC at u(z=0) = u(z=bin1))</span>
0757 <span class="comment">%     -There are cases where the bin corresponding with the bed actually</span>
0758 <span class="comment">%      contains flow data (i.e., not NaN or zero). For cases where the</span>
0759 <span class="comment">%      blanking distance DOES exists, I have imposed a BC of U=0 at the bed,</span>
0760 <span class="comment">%     -In cases where data goes all of the way to the bed, I have divided</span>
0761 <span class="comment">%      the last bin's velocity by 2, essentially imposing a U=0 at the</span>
0762 <span class="comment">%      boundary by extrapolating to the bottom of the bin.</span>
0763 <span class="comment">%     -This function still needs to be incorporated into the GUI.</span>
0764 
0765 <span class="comment">% V2 modifies the code for integration into the VMT GUI. Adds the Rozovskii output to</span>
0766 <span class="comment">% the V structure and computes the X components of the primary and secondary</span>
0767 <span class="comment">% velocities (in addition to cross stream Y components).</span>
0768 <span class="comment">% P.R. Jackson, USGS</span>
0769 
0770 
0771 <span class="comment">% Written by:</span>
0772 <span class="comment">% Frank L. Engel (fengel2@illinois.edu)</span>
0773 
0774 <span class="comment">% Last edited: 6/10/2010 FE</span>
0775 <span class="comment">% FE- I fixed an error in computing theta for data in quadrants 3 &amp; 4. The</span>
0776 <span class="comment">%     linear interpolation of velocity to the bed BC was creating errors</span>
0777 <span class="comment">%     in the computation of us, so I removed it. Also, I made a new</span>
0778 <span class="comment">%     variable which is the sum of all vs as an error check (it should</span>
0779 <span class="comment">%     always sum to zero)</span>
0780 
0781 disp(<span class="string">'Performing Rozovskii analysis...'</span>)
0782 bin_size = A(1,1).Sup.binSize_cm/100; <span class="comment">% in meters</span>
0783 
0784 <span class="keyword">for</span> i = 1:size(V.mcsMag,2)
0785     <span class="comment">% Finds closest bin to beam avg. depth (ie from V.mcsBed)</span>
0786     [min_difference(i), array_position(i)]<span class="keyword">...</span>
0787         = min(abs(V.mcsDepth(:,i) - V.mcsBed(i)));
0788     <span class="comment">%disp(['ap = ' num2str(array_position(i))])</span>
0789     <span class="comment">% Create a seperate version of the velocity data which can be modified,</span>
0790     <span class="comment">% preserving the VMT processing. Replaces all of the NaNs with u=0.</span>
0791     temp_u = V.u;
0792     temp_v = V.v;
0793     n = find(isnan(V.u));
0794     temp_u(n) = 0;
0795     temp_v(n) = 0;
0796     
0797     <span class="comment">% Compute Depth-averaged velocities and angles (using a difference</span>
0798     <span class="comment">% scheme)</span>
0799     <span class="keyword">for</span> j = 1:array_position(i)
0800         <span class="keyword">if</span> j == 1 <span class="comment">% Near water surface</span>
0801             <span class="comment">% Compute first bin by exprapolating velocity to the water</span>
0802             <span class="comment">% surface. WSE = 0. Imposes BC u(z=0) = u(z=bin1)</span>
0803             du_i(j,i) = temp_u(j,i)*(V.mcsDepth(j+1,i)-V.mcsDepth(j,i))<span class="keyword">...</span>
0804                       + temp_u(j,i)*(V.mcsDepth(j,i)-bin_size/2-0);
0805             dv_i(j,i) = temp_v(j,i)*(V.mcsDepth(j+1,i)-V.mcsDepth(j,i))<span class="keyword">...</span>
0806                       + temp_v(j,i)*(V.mcsDepth(j,i)-bin_size/2-0);
0807         <span class="keyword">elseif</span> j &lt; array_position(i) <span class="comment">% Inbetween</span>
0808             du_i(j,i) = temp_u(j,i)*(V.mcsDepth(j+1,i)-V.mcsDepth(j-1,i))/2;
0809             dv_i(j,i) = temp_v(j,i)*(V.mcsDepth(j+1,i)-V.mcsDepth(j-1,i))/2;
0810             
0811             <span class="comment">% Got rid of the linear interpolation to the bed- it was</span>
0812             <span class="comment">% messing up the Rozovskii secondary velocities (FE 6/10/2010)</span>
0813         <span class="keyword">elseif</span> j == array_position(i) <span class="comment">% Near bed</span>
0814             <span class="comment">%             k=0;</span>
0815             <span class="comment">%             % Search bins above the bed for the first bin containing flow</span>
0816             <span class="comment">%             % data.</span>
0817             <span class="comment">%</span>
0818             <span class="comment">% %             while temp_u(j-k,i) == 0</span>
0819             <span class="comment">% %                 if temp_u(j-k,i) == 0; k = k + 1; else end % find next good bin</span>
0820             <span class="comment">% %             end</span>
0821             <span class="comment">%</span>
0822             indx = find(temp_u(:,i) ~= 0);  <span class="comment">%Revision PRJ 9-1-09</span>
0823             <span class="keyword">if</span> isempty(indx)
0824                 du_i(:,i) = NaN;
0825                 dv_i(:,i) = NaN;
0826             <span class="keyword">else</span>
0827                 l = indx(end);
0828                 k = j - l;
0829                 <span class="comment">% Computes du from last good bin to the bed by linear</span>
0830                 <span class="comment">% interpolation. IMPOSES BC: u=0 at the bed</span>
0831 <span class="comment">%                 du_i(j-k+1,i) = (temp_u(j-k,i)-temp_u(j,i))/k...</span>
0832 <span class="comment">%                     *(V.mcsDepth(j,i)-V.mcsDepth(j-k,i))/k;</span>
0833 <span class="comment">%                 dv_i(j-k+1,i) = (temp_v(j-k,i)-temp_v(j,i))/k...</span>
0834 <span class="comment">%                     *(V.mcsDepth(j,i)-V.mcsDepth(j-k,i))/k;</span>
0835                 
0836                 <span class="comment">% Paints everything below last bin as NaN</span>
0837                 du_i(j-k+2:size(V.u,2),i) = NaN;
0838                 dv_i(j-k+2:size(V.u,2),i) = NaN;
0839             <span class="keyword">end</span>
0840         <span class="keyword">end</span>
0841     <span class="keyword">end</span>
0842     
0843     <span class="comment">% Depth averaged quantities</span>
0844     U(i) = nansum(du_i(:,i))/V.mcsDepth(array_position(i),i);
0845     V1(i) = nansum(dv_i(:,i))/V.mcsDepth(array_position(i),i);
0846     U_mag(i) = sqrt(U(i)^2+V1(i)^2); <span class="comment">% resultant vector</span>
0847     
0848     <span class="comment">% Angle of resultant vector from a perpendicular line along the</span>
0849     <span class="comment">% transect</span>
0850     phi(i) = atan(V1(i)/U(i));
0851     phi_deg(i) = rad2deg(phi(i));
0852     
0853     <span class="comment">% Compute Rozovskii variables at each bin</span>
0854     <span class="keyword">for</span> j = 1:array_position(i)
0855         u(j,i) = sqrt(V.u(j,i)^2+V.v(j,i)^2);
0856         <span class="keyword">if</span> (V.u(j,i) &lt; 0) &amp;&amp; (V.v(j,i) &lt; 0)
0857             theta(j,i) = atan(V.v(j,i)/V.u(j,i)) - pi();
0858         <span class="keyword">elseif</span> (V.u(j,i) &lt; 0) &amp;&amp; (V.v(j,i) &gt; 0)
0859             theta(j,i) = atan(V.v(j,i)/V.u(j,i)) + pi();
0860         <span class="keyword">else</span>
0861             theta(j,i) = atan(V.v(j,i)/V.u(j,i));
0862         <span class="keyword">end</span>
0863         theta_deg(j,i) = rad2deg(theta(j,i));
0864         up(j,i) = u(j,i)*cos(theta(j,i)-phi(i));
0865         us(j,i) = u(j,i)*sin(theta(j,i)-phi(i));
0866         upy(j,i) = up(j,i)*sin(phi(i));
0867         upx(j,i) = up(j,i)*cos(phi(i));
0868         usy(j,i) = us(j,i)*cos(phi(i));
0869         usx(j,i) = us(j,i)*sin(phi(i));
0870         depths(j,i) = V.mcsDepth(j,i);
0871         
0872         <span class="comment">% Compute d_us to check for zero secondary discharge constraint</span>
0873         <span class="keyword">if</span> j == 1 <span class="comment">% Near water surface</span>
0874             dus_i(j,i) = us(j,i)*(V.mcsDepth(j+1,i)-V.mcsDepth(j,i))<span class="keyword">...</span>
0875                 + us(j,i)*(V.mcsDepth(j,i)-bin_size/2-0);
0876         <span class="keyword">elseif</span> j &lt; array_position(i) <span class="comment">% Inbetween</span>
0877             dus_i(j,i) = us(j,i)*(V.mcsDepth(j+1,i)-V.mcsDepth(j-1,i))/2;
0878         <span class="keyword">end</span>
0879         <span class="comment">% Sum dus_i - this should be near zero for each ensemble</span>
0880         q_us(i) = nansum(dus_i(:,i));
0881     <span class="keyword">end</span>
0882     
0883     <span class="comment">% Resize variables to be the same as V structure array</span>
0884     indices = j+1:size(V.mcsMag,1);
0885     u(indices,i)         = NaN;
0886     theta(indices,i)     = NaN;
0887     theta_deg(indices,i) = NaN;
0888     up(indices,i)        = NaN;
0889     us(indices,i)        = NaN;
0890     upy(indices,i)       = NaN;
0891     usy(indices,i)       = NaN;
0892     upx(indices,i)       = NaN;
0893     usx(indices,i)       = NaN;
0894     depths(indices,i)    = NaN;
0895     dus_i(indices,i)     = NaN;
0896 <span class="comment">%     u(j+1:size(V.mcsMag,1),i) = NaN;</span>
0897 <span class="comment">%     theta(j+1:size(V.mcsMag,1),i) = NaN;</span>
0898 <span class="comment">%     theta_deg(j+1:size(V.mcsMag,1),i) = NaN;</span>
0899 <span class="comment">%     up(j+1:size(V.mcsMag,1),i) = NaN;</span>
0900 <span class="comment">%     us(j+1:size(V.mcsMag,1),i) = NaN;</span>
0901 <span class="comment">%     upy(j+1:size(V.mcsMag,1),i) = NaN;</span>
0902 <span class="comment">%     usy(j+1:size(V.mcsMag,1),i) = NaN;</span>
0903 <span class="comment">%     upx(j+1:size(V.mcsMag,1),i) = NaN;</span>
0904 <span class="comment">%     usx(j+1:size(V.mcsMag,1),i) = NaN;</span>
0905 <span class="comment">%     depths(j+1:size(V.mcsMag,1),i) = NaN;</span>
0906 <span class="comment">%     dus_i(j+1:size(V.mcsMag,1),i) = NaN;</span>
0907 <span class="keyword">end</span>
0908 
0909 <span class="comment">% Display error message if rozovskii computation of q_us doesn't sum to</span>
0910 <span class="comment">% zero</span>
0911 <span class="keyword">if</span> q_us &gt; 1e-4
0912     disp(<span class="string">'Warning: Rozovskii secondary velocities not satisfying continuity!'</span>)
0913 <span class="keyword">else</span>
0914     disp(<span class="string">'Computation successfull: Rozovskii secondary velocities satisfy continuity'</span>)
0915 <span class="keyword">end</span>
0916 
0917 <span class="comment">% Rotate local velocity vectors into global coordinate system by</span>
0918 <span class="comment">% determining the angle of the transect using endpoint locations. The</span>
0919 <span class="comment">% function &quot;vrotation&quot; is a standard rotation matrix</span>
0920 XStheta = atan((V.mcsY(1,end)-V.mcsY(1,1))/(V.mcsX(1,end)-V.mcsX(1,1)));
0921 XSalpha = XStheta - pi/2;
0922 [ux, uy, uz] = vrotation(V.u,V.v,V.w,XSalpha);
0923 
0924 <span class="comment">% Put results into the V structure</span>
0925 
0926 V.Roz.U         = U;
0927 V.Roz.V         = V1;
0928 V.Roz.U_mag     = U_mag;
0929 V.Roz.phi       = phi;
0930 V.Roz.phi_deg   = phi_deg;
0931 V.Roz.u         = V.u;
0932 V.Roz.v         = V.v;
0933 V.Roz.u_mag     = u;
0934 V.Roz.depth     = depths;
0935 V.Roz.theta     = theta;
0936 V.Roz.theta_deg = theta_deg;
0937 V.Roz.up        = up;
0938 V.Roz.us        = us;
0939 V.Roz.upy       = upy;
0940 V.Roz.usy       = usy;
0941 V.Roz.upx       = upx;
0942 V.Roz.usx       = usx;
0943 V.Roz.ux        = ux;
0944 V.Roz.uy        = uy;
0945 V.Roz.uz        = uz;
0946 V.Roz.alpha     = XSalpha;
0947 
0948 <span class="comment">% Rozovskii = struct('Ux', {Ux}, 'Uy', {Uy}, 'U', {U}, 'phi', {phi},...</span>
0949 <span class="comment">%     'phi_deg', {phi_deg},'ux', {V.u}, 'uy', {V.v}, 'u', {u}, 'depth', {depths},...</span>
0950 <span class="comment">%     'theta', {theta}, 'theta_deg', {theta_deg}, 'up', {up}, 'us', {us},...</span>
0951 <span class="comment">%     'upy', {upy}, 'usy', {usy});</span>
0952 
0953 <span class="comment">% Name of output file needs to be modified to take handle args from GUI</span>
0954 <span class="comment">% Save variable into the VMTProcFile folder</span>
0955 <span class="comment">% outfolder = ['VMTProcFiles\'];</span>
0956 <span class="comment">% outfile=['Rozovskii'];</span>
0957 <span class="comment">% filename = [outfolder outfile];</span>
0958 <span class="comment">% save(filename, 'Rozovskii');</span>
0959 
0960 disp(<span class="string">'Rozovskii analysis complete. Added .Roz structure to V data structure.'</span>)
0961 <span class="comment">% directory = pwd;</span>
0962 <span class="comment">% fileloc = [directory '\' filename '.mat'];</span>
0963 <span class="comment">% disp(fileloc)</span>
0964</pre></div>
<hr><address>Generated on Thu 21-Mar-2013 09:32:01 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" target="_parent">m2html</a></strong> &copy; 2005</address>
</body>
</html>